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Abstract 

In Multi-Input Multi-Output (MIMO) systems, Maximum-Likelihood (ML) decoding is 
equivalent to finding the closest lattice point in an A^-dimensional complex space. In general, 
this problem is known to be NP hard. In this paper, we propose a quasi-maximum likelihood 
algorithm based on Semi-Definite Programming (SDP). We introduce several SDP relaxation 
models for MIMO systems, with increasing complexity. We use interior-point methods for solv- 
ing the models and obtain a near-ML performance with polynomial computational complexity. 
Lattice basis reduction is applied to further reduce the computational complexity of solving 
these models. The proposed relaxation models are also used for soft output decoding in MIMO 
systems. 
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I. Introduction 

Recently, there has been a considerable interest in Multi-Input Multi-Output (MIMO) 
antenna systems due to achieving a very high capacity as compared to single- antenna 
systems [4], [5]. In MIMO systems, a vector is transmitted by the transmit antennas. In 
the receiver, a corrupted version of this vector affected by the channel noise and fading is 
received. Decoding concerns the operation of recovering the transmitted vector from the 
received signal. This problem is usually expressed in terms of "lattice decoding" which 
is known to be NP-hard. 

To overcome the complexity issue, a variety of sub-optimum polynomial time algo- 
rithms are suggested in the literature for lattice decoding. However, unfortunately, these 
algorithms usually result in a noticeable degradation in the performance. Examples of 
such polynomial time algorithms include: Zero Forcing Detector (ZFD) [6], [7], Minimum 
Mean Squared Error Detector (MMSED) [8], [9], Decision Feedback Detector (DFD) 
and Vertical Bell Laboratories Layered Space-Time Nulling and Cancellation Detector 
(VBLAST Detector) [10], [11]. 

Lattice basis reduction has been applied as a pre-processing step in sub-optimum de- 
coding algorithms to reduce the complexity and achieve a better performance. Minkowski 
reduction [12], Korkin-Zolotarev reduction [13] and LLL reduction [14] have been suc- 
cessfully used for this purpose in [15-20]. 

In the last decade. Sphere Decoder (SdJi] is introduced as a Maximum Likelihood 
(ML) decoding method for MIMO systems with near-optimal performance [23]. In the 
SD method, the lattice points inside a hyper-sphere are generated and the closest lattice 
point to the received signal is determined. The average complexity of this algorithm 
is shown to be polynomial time (almost cubic) over certain ranges of rate. Signal to 
Noise Ratio (SNR) and dimension, while the worst case complexity is still exponential 
[15], [24], [25]. However, recently, it has been shown that it is a misconception that the 
expected number of operations in SD asymptotically grows as a polynomial function of 

'This technique is introduced in the mathematical literature several years ago [21], [22]. 
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the problem size [26] (reference [26] derives an exponential lower bound on the average 
complexity of SD). 

In [27], a quasi-maximum likelihood method for lattice decoding is introduced. 
Each signal constellation is expressed by its binary representation and the decoding is 
transformed into a quadratic minimization problem [27]. Then, the resulting problem is 
solved using a relaxation for rank-one matrices in Semi-Definite Programming (SDP) 
context. It is shown that this method has a near optimum performance and a polynomial 
time worst case complexity. However, the method proposed in [27] is limited to scenarios 
that the constellation points are expressed as a linear combination of bit labels. A 
typical example is the case of natural labeling in conjunction with PSK constellation. 
This restriction is removed in this work. Therefore, we are able to handle any form of 
constellation for an arbitrary labeling of points, for example PAM constellation with Gray 
labeling^. Another quasi-maximum likelihood decoding method is introduced in [29] for 
larger PSK constellations with near ML performance and low complexity. 

Recently, we became aware of another quasi-maximum likelihood decoding method 
[30] for the MIMO systems employing 16-QAM. They replace any finite constellation 
by a polynomial constraint, e.g. if x G {a, b, c}, then (x — a)(x — h){x — c) = 0. Then, 
by introducing some slack variables, the constraints are expressed in terms of quadratic 
polynomials. Finally, the SDP relaxation is resulted by dropping the rank-one constraint. 
The work in [30] , in its current form, is restricted to MIMO systems employing 16-QAM. 
However, it can be generalized for larger constellations at the cost of defining more slack 
variables, increasing the complexity, and significantly decreasing the performance^. 

In this work, we develop an efficient approximate ML decoder for MIMO systems 
based on SDP. In the proposed method, the transmitted vector is expanded as a linear 
combination (with zero-one coefficients) of all the possible constellation points in each 
dimension. Using this formulation, the distance minimization in Euclidean space is ex- 

^It is shown that Gray labeling, among all possible constellation labeling methods, offers the lowest possible average 
probability of bit errors [28]. 
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pressed in terms of a binary quadratic minimization problem. The minimization of this 
problem is over the set of all binary rank-one matrices with column sums equal to one. In 
order to solve this minimization problem, we present two relaxation models (Model HI 
and IV), providing a trade-off between the computational complexity and the performance 
(both models can be solved with polynomial-time complexity). Two additional relaxation 
models (Model I and II) are presented as intermediate steps in the derivations of Model 
m and IV. 

Model I: A preliminary SDP relaxation of the minimization problem is obtained 
by removing the rank-one constraint in the problem and using Lagrangian duality [31]. 
This relaxation has many redundant constraints and no strict interior point in the feasible 
set (there are numerical difficulties in computing the solution for a problem without an 
interior point). 

Model II: To overcome this drawback, the feasible set is projected onto a face of 
the semi-definite cone. Then, based on the identified redundant constraints, another form 
of the relaxation is obtained, which can be solved using interior-point methods. 

Model III: The relaxation Model II results in a weak lower bound. To strengthen this 
relaxation model, the structure of the feasible set is investigated. An interesting property 
of the feasible set imposes a zero pattern for the solution. Adding this pattern as an extra 
constraint to the previous relaxation model results in a stronger model. 

Model IV: Finally, the strongest relaxation model in this work is introduced by 
adding some additional non-negativity constraints. The number of non-negativity con- 
straints can be adjusted to provide a trade-off between the performance and complexity 
of the resulting method. 

Simulation results show that the performance of the last model is near optimal for 
M-ary QAM or PSK constellation (with an arbitrary binary labeling, say Gray labeling). 
Therefore, the decoding algorithm built on this model has a near-ML performance with 
polynomial computational complexity. 

The proposed models result in a solution that is not necessarily a binary rank-one 
matrix. This solution is changed to a binary rank-one matrix through a randomization 
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algorithm. We modify the conventional randomization algorithms to adopt to our prob- 
lem. Also, a new randomization procedure is introduced which finds the optimal binary 
rank-one solution in a smaller number of iterations than the conventional ones. Finally, 
we discuss using a lattice basis reduction method to further reduce the computational 
complexity of the proposed relaxation models. The extension of the decoding technique 
for soft output decoding is also investigated. 

Following notations are used in the sequel. The space of K x N (resp. N x N) 
real matrices is denoted by Mkxn (resp. A4n), and the space of N x N symmetric 
matrices is denoted by Sn- The indexing of the matrix elements start from one, unless 
in cases that an extra row (or column) is inserted in the first row (or column) of a matrix 
in which case indexing starts from zero. For a K x N matrix X e Mrxn the {i,j)th 
element is represented by Xij, where I < i < K, 1 < j < N, i.e. X = [xij]. We use 
tracc(A) to denote the trace of a square matrix A. The space of symmetric matrices is 
considered with the trace inner product ( A, B) = tracc(AB). For A, B G Sn, A ^ 
(resp. A 0) denotes positive semi-definiteness (resp. positive definiteness), and A ^ B 
denotes A — B >: 0. For two matrices A, B e M.n, A > B, (A > B) means > by, 
{ttij > bij) for all The Kronecker product of two matrices A and B is denoted by 
A B (for definition see [32]). 

For X e M.KxN, vec(X) denotes the vector in (real XA^-dimensional space) 
that is formed from the columns of the matrix X. The following identity relates the 
Kronecker product with vec( ) operator, see e.g. [32], 

vec( ACB) = (B^ ® A) vec(C) . (1) 

For X G TWat, diag(X) is a vector of the diagonal elements of X. We use e^r G (resp. 
Oat G R^) to denote the x 1 vector of all ones (resp. all zeros), ^k-kn & Mrxn 
to denote the matrix of all ones, and Ijv to denote the N x N Identity matrix. For 
X e M-KxN, the notation X(l : i,l : j), i < k and j < n denotes the sub-matrix of X 
containing the first i rows and the first j columns. 

The rest of the paper is organized as follows. The problem formulation is introduced 
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in Section II. Section HI is devoted to the semi-definite solution of this problem. In 
Section IV, randomization procedures for finding rank-one solutions are presented. Section 
V introduces a method based on lattice basis reduction to reduce the computational 
complexity of the proposed relaxation models. In Section VI, the soft decoding methods 
based on the proposed models are investigated. Simulation results are presented in Section 
VII. Finally, Section VIII concludes the paper. 

II. Problem Formulation 
A MIMO system with N transmit antennas and M receive antennas is modelled as 



^'^^m + fl, (2) 



ME, 



where H 



hij 



is the M X N channel matrix composed of independent, identically 
distributed complex Gaussian random elements with zero mean and unit variance, n 
is an M X 1 complex additive white Gaussian noise vector with zero mean and unit 
variance, and x is an iV x 1 data vector whose components are selected from a complex 
set (S = {si, S2, ■ ■ ■ , sj^} with an average energy of Eg^^. The parameter SNR in ^ is 
the SNR per receive antenna. 

Noting G iS, for i = 1, ■ ■ ■ , iV, we have 

Xi = Ui{l)si + Ui{2)s2 -\ ^Ui{K)sK, (3) 

where 

,(j)e{0,l} and ^M,(j) = l, V2 = l,---,iV. (4) 



Let u 



K 

1 T 



^ii(l) •■■ ui{K) ■■■ UNiX) ■■■ MK) 



and N = N. Using the equa- 



tions in ([3]) and dH), the transmitted vector is expressed as 

± = Su, (5) 

where S = Iat (g) [si, ■ ■ ■ , sk] is an x NK matrix of coefficients, and u is an NK x 1 
binary vector such that Au = ejv, where A = Ijv ® e]^. This constraint states that 
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among each K components of the binary vector u, i.e. ■ ■ ■ ,Ui{K), there is only 

one element equal to "1" and the rest are zero. 

To avoid using complex matrices, the system model (|2l) is represented by real 
matrices in 

-J(H) 9^(H 



^(y) 
3f(y) 



+ 



MEs 

J(n)_ 
Hx + n 



3{±) 



(6) 



where 9^(.) and 3(.) denote the real and imaginary parts of a matrix, respectively, y is 
the received vector, and x is the input vector. 



Let S denotes the real matrix 



; therefore. 



y = HSu + n 



(7) 



expresses the MIMO system model by real matrices and the input binary data vector, u. 

Consider the case that different components of x in Q, corresponding to the two- 
dimensional sub-constellations, are equal to the cartesian product of their underlying 
one-dimensional sub-constellations, e.g. QAM signalling. In this case, the components 
of X in (l6l) belong to the set S = {si, ■ ■ ■ , sk} with real elements, i.e. 



Xi = Ui{l)si + Ui{2)s2 H h Ui{K)sK, 



(8) 



where only one of the Ui{j) is 1 and the rest are zero. 

Let u = [ni(l) ■■■ui{K)--- un{1) ■ ■ ■UN{K)f, N = 2N, and S = lAr®[si, ■ ■ ■ ,sk] 
Then, the equation for the components of x in ([8]) reduces to x = Su and the relationship 
for the MIMO system model is given in ([7]). 

At the receiver, the Maximum-Likelihood (ML) decoding rule is given by 



X = argmin ||y — Hx| 



(9) 
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where x is the most likely input vector and y is the received vector. Noting x = Su, 
this problem is equivalent to 

min ||y — HSu|p = 

Au=e]v 

min u^S^H^HSu - 2y'^HSu, (10) 

Au=ejv 

where u is a binary vector. 

Let Q = S^H'^HS and c = — S^H^y. Therefore, this problem is formulated as 

min u^Qu + 2c^u 

S.t. Au = Gat 



where n = NK. The formulation (fTTl) is a quadratic minimization problem with binary 
variables [31]. Recent studies on solving binary quadratic minimization problems such as 
Graph Partitioning [33] and Quadratic Assignment Problem [34], [35], show that semi- 
definite programming is a very promising approach to provide tight relaxations for such 
problems. In the following, we derive several SDP relaxation models for the minimization 
problem in (fTT)) . Appendix H] provides the mathematical framework for these models using 
the Lagrangian Duality [31]. 



Consider the minimization problem in (fTT]) . Since u is a binary vector, the objective 
function is expressed as 



u G {0,1} 



(11) 



in. Semi-Definite Programming Solution 




(12) 
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where £q := 

Let 8k-> 
to one, i.e. 








c 


Q 


denote f 



le set of all binary matrices in M.kxn with column sums equal 



£KxN = {y^^MKxN ■ e^X = e^,Xjj G {0, l},Vi,j} . 



(13) 



Since the constraints Au = e^r, Ui G {0, 1}^^' in (ITTI) and u = vec(U), U G Ekxn are 
equivalent, the minimization problem (fTTI) can be written as 



min trace C 



Q 



1 


T 

U 


u 


T 

UU 



(14) 



s.t. u = vec(U), U G £kxN- 
To derive the first semi-definite relaxation model, a direct approach based on the 
well known lifting process [36] is selected. In accordance to (fT4l) . for any U G Srxn, 
u = vec(U), the feasible points of (fT4l) are expressed by 



1 
u 



1 



1 



u 



u 



UU 



(15) 



The matrix Yu is a rank-one and positive semi-definite matrix. Also, we have 

diag(Yj = Y;^^_^ = Y„ 



where Y^q . (resp. Yu g) denotes the first row (resp. the first column]^ of Yu (Note that 
u is a binary vector, and consequently, diag(uu^) = u). 

In order to obtain a tractable SDP relaxation of (fT4l) . we remove the rank-one 
restriction from the feasible set. In fact, the feasible set is approximated by another 
larger set JF, defined as 



:= conv {Yu : u = vec(U), U G Srxn} 

''Matrix Yu is indexed from zero. 



(16) 
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where conv(.) denotes the convex hull of a set. This results in our first relaxation model 
(Model I) for the original problem given in (fTTI) : 

min trace^oY 

^ (17) 

s.t. Y G 

It is clear that the matrices 

Yu for u = vec(U), U G Skxn 

are the feasible points of JF. Moreover, since these points are rank-one matrices, they are 
contained in the set of extreme points of JF, see e.g. [37]. In other words, if the matrix 

_ 1 
Y is restricted to be rank-one in (flTl) . i.e. Y = 

the optimal solution of (flTl) provides the optimal 



u 



[l u^] , for some u G M", then 



solution of (fTTI) . 

The SDP relaxation problem (fTTI) is not solvable in polynomial time and JF has no 
interior points. Therefore our goal is to approximate the set JF by a larger set containing 
JF. In the following, we show that JF actually lies in a smaller dimensional subspace. We 
will further see that relative to this subspace, JF will have interior points. 



A. Geometry of the Relaxation 

In order to approximate the feasible set JF for solving the problem, we elaborate 
more on the geometrical structure of this set. First, we prove the following lemma on 
the representation of matrices having sum of the elements in each column equal to one. 

Lemma 1: Let 



V 



Kx{K-l) 



G M 



Kx{K-l) 



and 



KxN 



K L 



KxN 



Kx{K-l)^(K-l)xN 



(18) 



(19) 



A matrix X G A4kxn with the property that the summation of its elements in each 
column is equal to one, i.e. e]^X = e^, can be written as 



X — Fkxn + Vi^x(A'-i)Z, 



(20) 
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where Z = X(l : {K-l),l : A^). 

Proof: see Appendix Ull ■ 
Corollary 1: VX G £kxN, 3Z G M(k~i)xN, G {0,1} s.t. X = 

where Z = X(l : {K : A^). 

Using Lemma [H the following theorem can be proved which provides the structure of 
the elements in the set JF. 
Theorem 2: Let 



V 



1 


0^ 

^N(K-l) 


-^{^NK — (lAf ® ^ Kx(K-^1))^(K-1)n) 





(21) 



where V G M.[nk+i)x({k-i)n+i)- For any Y G JF, there exists a symmetric matrix R of 
order N{K - 1) + 1, indexed from to N{K - 1), such that 

Y = VRV^, R ^ 0, and roo = 1, r,, = r^i, (22) 

Also, if Y is an extreme point of JF, then r^j G {0, 1}, otherwise Vij G [0, 1] for i,j G 
{0,...,iV(i^-l)}. 

Proof: see Appendix |Ill ■ 
Using Theorem [21 we can show that the set JF^ contains JF: 

jFj, = |Y G iSata'+i : 3 R G S(^K~i)N+i, R ^ 0, 

/2oo = 1, Y = VRV^, diag(Y) = Yq,} . (23) 

Therefore, the feasible set in (fTT] ) is approximated by JF,,. This results in our second 
relaxation model (Model II) of the original problem given in (fTTI) : 

min trace (V^£qV)R 
s.t. diag(VRV^) = (1, (VRv'^)o,i;„)^ 

R ^ 0. (24) 

Note that the matrices Yu are contained in the set of extreme points of JF. We need 
only consider faces of JF which contain all of these extreme points. Therefore, we are 
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only looking for the minimal face, which is the intersection of all these faces. We will 
show later that the SDP relaxation (|24l) is the projection of the SDP relaxation (fTT]) onto 
the minimal face of JF. 

Solving the relaxation model in (flTI) over JF results in the optimal solution of the 
original problem in (fT4l) . but this problem is NP-hard. Solving the relaxation model in 
(|24|) over JF^ results in a weaker bound for the optimal solution. In order to improve this 
bound, the relaxation is strengthen by adding an interesting property of the matrix Yu. 
This results in the next relaxation model. 

B. Tightening the Relaxation by Gangster Operator 

The feasible set of the minimization problem (l24l) is convex. It contains the set 
of matrices of the form Yu corresponding to different vectors u. However, the SDP 
relaxations may contain many points that are not in the affine hull of these Yu. In the 
following, we extract a condition which is implicit in the matrix Yu and explicitly add it 
to the relaxation model (l24l) . Subsequently, some redundant constraint are removed and 
this results in an improved relaxation (relaxation Model III). 

Theorem 3: Let U denote the set of all binary vectors u = vec(U), U G Srxn- 
Define the barycenter point, Y, as the arithmetic mean of all the feasible points in the 
minimization problem (fT4l) : therefore. 



i) Y has (a) the value of 1 as its (0, 0) element, (b) blocks of dimension K x K 
on its diagonal which are diagonal matrices with elements 1/K, and (c) the first 
row and first column equal to the vector of its diagonal elements. The rest of the 




(25) 



Then: 
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matrix is composed of K x K blocks with all elements equal to l/K"^ 



1 




1 T 
— P 






1 T 























1 T 
— P 



+ 



On 







(KIk - Ex) 



(26) 



ii) rank(Y) = N{K - 1) + 1; 



iii) The NK + 1 eigenvalues of Y are given in the vector 

iv) The null space of Y can be expressed by A/'(Y) = |m : m G 7^(T^)} , where the 
constraint matrix T is the following N x (NK + 1) matrix 



A 



v) the range of Y can be expressed by the columns of the (NK + 1) x (N(K — 1) + 1) 
matrix V. Furthermore, TV = 0. 

Proof: see Appendix UIl ■ 
Remark 2: The faces of the positive semi-definite cone are characterized by the 
null space of the points in their relative interior. The minimal face of the SDP problem 
contains matrices Yu and can be expressed as y'S]\i(^K-i)+i^'^ ■ Thus, the SDP relaxation 
(l24l) is a projected relaxation onto the minimal face of the feasible set JF. 

Theorem [3] suggests a zero pattern for the elements of JF. We use a Gangster 
Operator [34] to represent these constraints more efficiently. Let J be a set of indices, 
then this operator is defined as 
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Yij if or e J 
otherwise. 



(27) 



Considering the barycenter point, we have ^j(Y) = for 

J={{i,j): i = K{j>-l) + q, j = K{p-l)+r, 
q<r,q,r e {!,■■■ ,K},pe {!,■■■ , N}} . 



(28) 



Since Y is a convex combination of all matrices in U with entries either or 1; hence, 
from (|28l) . we have ^j(Yu) = 0. Also, all the points from the feasible set JF are the 
convex combination of Yu. Therefore, 

GjiY) = 0, VY G (29) 

The feasible set of the projected SDP in (l24l) is tightened by adding the constraints 
Qj(y) = 0. By combining these constraints and (|24|) . we note that there are some 
redundant constraints that can be removed to enhance the relaxation model. This is 
expressed in the following lemma. 

Lemma 4: Let R be an arbitrary {N{K — 1) + 1) x {N{K — 1) + 1) symmetric 
matrix with 





^00 


R-oi 


RoAf 


R = 


R-io 


Rii 


RlAT 










Rati 


Unn 



(30) 



where tqq is a scalar, Rjo, for z = 1, ■ ■ ■ , are (K — 1) x 1 vectors and Rjj, for 
i,j = !,■■■ ,N, are {K-l)x (K-l) blocks of R. Theorem [2] states that Y = VRV^. 
We can partition Y as 





Yol • 


Yqat 


Yio 


Yn • 


YlAT 


Yatq 


Yati 


Ytvat 



(31) 
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where ?/oo is a scalar, Yjo, for i = 1, ■ ■ ■ , N are Kxl vectors and Y^j, for i,j = 1, 
are K x K blocks of Y. Then, 

1) yoQ = roo and Yo^ei^ = roo, for z = 1, ■ 



N, 



N. 



2) Y, 



e^Yij for i,j = 1, 



AT. 



Proof: Noting TY = (see Theorem [3l 4), the proof follows. ■ 
If the Gangster operator is applied to ((24l) . it results in the following redundant constraint 



diag(VRV ) = (1,(VRV )o,l:„)^ 

Note that using Lemma HI Yqj = e^Yjj for j = 1, ■ ■ ■ , and the off-diagonal entries 
of each Yjj are zero. Therefore, by defining a new set J = JU{0,0} and eliminating 
the redundant constraints, we obtain a new SDP relaxation model (Model III): 

min tmce(V^ CqV)R 
s.t. 6?j(VRV^) = Eoo 

R ^ 0, (32) 

where R is an {N{K-l) + l)x{N{K-l) + l) matrix and Eqo is an {NK+l)x{NK+l) 
all zero matrix except for a single element equal to 1 in its (0, 0)th entry. With this new 
index set J, we are able to remove all the redundant constraints while maintaining the 
SDP relaxation. The relaxation model in (|32l ) corresponds to a tighter lower bound and 
has an interior point in its feasible set as shown in the following theorem. 
Theorem 5: The {N{K - 1) + 1) x {N{K - 1) + 1) matrix 



R 



K*^N{K-1) 



E 



K-l) 



(33) 



is a strictly interior point of the feasible set for the relaxation problem (1321) . 

Proof: The matrix R is positive definite. The rest of the proof follows by showing 
VRV^ = Y. ■ 

The relaxation in (l32l) is further tightened by considering the non-negativity con- 
straints [35]. All the elements of the matrix Y which are not covered by the Gangster 
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operator are greater than or equal to zero. These inequalities can be added to the set of 
constraints in (|32l) . resulting in a stronger relaxation model (Model IV): 

min trace(V^i:QV)R 
s.t. gjCVKV^) = Eoo 
^j(VRV^) > 

R ^ 0, (34) 

where the set J indicates those indices which are not covered by J. 

Note that this model is considerably stronger than model (|32|) because non-negativity 
constraints are also imposed in the model. The advantage of this formulation is that the 
number of inequalities can be adjusted to provide a trade-off between the strength of the 
bounds and the complexity of the problem. The larger number of the constraints in the 
model is, the better it approximates the optimization problem (fT4l) (with an increase in 
the complexity). 

The most common methods for solving SDP problems of moderate sizes (with 
dimensions on the order of hundreds) are IPMs, whose computational complexities are 
polynomial, see e.g. [38]. There are a large number of IPM -based solvers to handle SDP 
problems, e.g., DSDP [39], SeDuMi [40], SDPA [41], etc. In our numerical experiments, 
we use DSDP and SDPA for solving (l32l) . and SeDuMi is implemented for solving (|34|) . 
Note that adding the non-negativity constraints increases the computational complexity 
of the model. Since the problem sizes of our interest are moderate, the complexity of 
solving (|34l ) with IPM solvers is tractable. 

IV. Randomization Method 

Solving the SDP relaxation models (|32l) and (l34l) results in a matrix R. This matrix 
is transformed to Y using Y = VRV^, whose elements are between and 1. This 
matrix has to be converted to a binary rank-one solution of (fT4)) . i.e. Yu, or equivalently, 
a binary vector u as a solution for (fTTI) . 
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For any feasible point of (fT4l) . i.e. Yu, the first row, the first column, and the vector 
of the diagonal elements of this symmetric matrix are equal to a binary solution for (fTTI) . 
For any matrix Y resulting from the relaxation problems (l32l) or (|34l) . its first row, its 
first column, and the vector of its diagonal elements are equal. Therefore, the vector u is 
approximated by rounding off the elements of the first column of the matrix Y. However, 
this transformation results in a loose upper bound on the performance. In order to improve 
the performance, Y is transformed to a binary rank-one matrix through a randomization 
procedure. An intuitive explanation of the randomization procedure is presented in [27]. 
We present two randomization algorithms in the Appendix Hill 

V. Complexity Reduction Using Lattice Basis Reduction 

Lattice structures have been used frequently in different communication applications 
such as quantization or MIMO decoding. A real lattice A is a discrete set of M- 
dimensional vectors in the real Euclidean Af-space, M*"^, that forms a group under the 
ordinary vector addition. Every lattice A is generated by the integer linear combinations 
of a set of linearly independent vectors {bi, ■ ■ ■ , bjv}, where bj G A, and the integer 
N,N < M, is called the dimension of the latticqj. The set of vectors {bi, ■ ■ ■ ,h]\f} is 
called a basis of A, and the N x M matrix B = [bi, • ■ ■ ,bjv]^ which has the basis 
vectors as its rows is called the basis matrix (or generator matrix) of A. 

The basis for representing a lattice is not unique. Usually a basis consisting of 
relatively short and nearly orthogonal vectors is desirable. The procedure of finding such 
a basis for a lattice is called Lattice Basis Reduction. Several distinct notions of reduction 
have been studied, including Lenstra-Lenstra and Lovasz (LLL) reduced basis [14], which 
can be computed in polynomial time. 

An initial solution for the lattice decoding problem can be computed using one of the 
simple sub-optimal algorithms such as ZED or channel inversion, e.g. s' = [H^^y]. If the 
channel is not ill-conditioned, i.e. the columns of the channel matrix are nearly orthogonal 
and short, it is most likely that the ML solution of the lattice decoding problem is around 

'without loss of generality, we assume that A*' — M. 
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s'. Therefore, using a reduced basis for the lattice, each Xi in ^ can be expressed by a 
few points in S around s[, not all the points in S. In general, this results in a sub-optimal 
algorithm. However, for the special case of a MEMO system with two antennas (with real 
coefficients), it has been shown that by using the LLL approximation and considering 
two points per dimension we achieve the ML decoding performance [42]. 

Let L = HQ be the LLL reduced basis for the channel matrix H, where Q is a 
unimodular matrix. The MIMO system model in ^ can be written as 

y = LQ^x + n. (35) 

Consider the QAM signaling. Without loss of generality, we can assume coordinates of 
X are in the integer grid. Since Q is a unimodular matrix, the coordinates of a new 
variable defined as x' = Q^^x are also in the integer grid. Therefore, the system in 
(l35l) is modelled by y = Lx' + n. Note that by multiplying x by Q^^ the constellation 
boundary will change. However, it is shown that in the lattice decoding problem with 
finite constellations the best approach is to ignore the boundary and compute the solution 
[43]. If the solution is outside the region, it is considered as an error. This change of 
boundary will result in some performance degradation. The performance degradation for 
some scenarios are depicted in Fig. [5] and Fig. [6l 

In order to implement the proposed method using LLL basis reduction, each com- 
ponent of x' is expressed by a linear combination (with zero-one coefficients) of L 
(usually much smaller than K) integers around s'^, where s' = [L^ V]- Then, the proposed 
algorithm can be applied to this new model. Due to the change of constellation boundary, 
there is a degradation in the performance. However, the complexity reduction is large. The 
trade-off between performance degradation and complexity reduction can be controlled 
by the choice of L (see simulation results). The reduction in the complexity is more 
pronounced for larger constellations. Note that the dimension of the semi-definite matrix 
Y is * (K — 1) + 1. Therefore, the LLL reduction decreases the dimension of the 
matrix Y to * (L — 1) + 1 (where usually L <^ K), and consequently, decreases the 
computational complexity of the proposed algorithm. The performance of this method is 
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shown in the simulation results. 



VI. Extension for Soft Decoding 

In this section, we extend our proposed SDP relaxation decoding method for soft 
decoding in MIMO systems. The SDP soft decoder is derived as an efficient solution 
of the max- approximated soft ML decoder. The complexity of this method is much less 
than that in the soft ML decoder. Moreover, the performance of the proposed method is 
comparable with that in the ML one. Also, the proposed method can be applied to any 
arbitrary constellation and labelling method, say Grey labeling. 

In the MIMO system defined in any transmit data x is represented by A*";, = 
log2 K bits (x = map(b), where b is the corresponding binary input). Given a received 
vector y, the soft decoder returns the soft information about the likelihood of bj = 
or 1, j = 1, ■ ■ ■ , A^A*";,. The likelihoods are calculated by Log-Likelihood Ratios (LLR) 
in a Maximum A Posteriori (MAP) decoder by 



Define 



CAW = ^o,gi0. (37) 



It is shown that the LLR values are formulated by [44] 

EbGB,,i P(y|b). exp (^ibJ].L^,[fc] 



C{bj\y) = log 



EbeB,,oP(y|b)-exp (|bJ].LA,[fc] 



+ 'CAib.ly), (38) 

where h[k] denotes the sub- vector of b obtained by omitting its A;th element b^, i-iA,[k] 
denotes the vector of all Ca values, also omitting bk, and Bfc,i (resp. Bfc^o) denotes the set 
of all input vectors, b, such that bk = I (resp. bk = 0). Note that there is an isomorphism 
between B^^i (resp. B^^o) and X^^i (resp. Xk^), where X^^i (resp. X^^q) denotes the set 
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of all corresponding constellation symbols, i = {x : x = map(b),b G M^^i} (resp. 
Xfc^o = {x : X = map(b), b G Bfc,o})- 

As shown in [44], the computation of the LLR values in (l38l) requires computing 
the likelihood function p(y|b), i.e. 

p(y|x = map(b)) = hl^- '1^:^^ , (39) 

where cx^ 



SNR- 

By having the likelihood functions, these LLR values are approximated efficiently 
using the Max-log approximation [44] 

^E{bk\y) ^ + i maxl-^ || y-Hx f +bj^].L^,[fe]l 

Without loss of generality, we assume that all components, Xi, of an input vector are 
equiprobableO; therefore, the second term in each maximization in ( |40l ) will be removed. 
Hence, computing the LLR values requires to solve problems of the form 

min II y- Hx f , (41) 

where A; = 1, ■ ■ ■ , A^A*";, and C = or 1. Note that, as mentioned in [46], only NNi, + 1 
problems among 2NNb problems of the form (|4TI) are considered. 

The Quasi-Maximum likelihood decoding method proposed in this paper can be ap- 
plied to the problem (|4TI) . However, X^^ must be defined in implementing the algorithm. 
This set includes all the input vectors, x G , such that bk = C- Assigning or 1 to 
one of the bits in b removes half of the points in . In other words, when = one 
of the components of the input vector x, say Xp, can only select half of the points in the 
set S, say |sp^, ■ ■ ■ , Sp^^ |. Therefore, the pth component of x is represented by 

Xp = Up{l)sp^^ ^-Mp(^)sp^. (42) 

*In order to consider the effects of non-equiprobable symbols, both approaches presented in [45] and [46] can be 
apphed. 
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As a result, we have the same matrix expression x = Su as ([8]), except that the 
length of the vector uis (N — 1) * K + ^ and in the pth row of the matrix S, we have 
y elements |sp^, • • ■ , Sp^^ |, instead of K elements {si, ■ ■ ■ ,sk}- Now, the proposed 
method can be applied to the new equation based on the new matrix S and u. 

VII. Simulation Results 

A. Performance Analysis 

We simulate the two proposed models (|32|) and (|34|) for decoding in MIMO systems 
with QAM and PSK constellations. Fig. [T] demonstrates that the proposed quasi-ML 
method using model (|32l) and the randomization procedure achieves near ML perfor- 
mance in an un-coded 2x2 MIMO system with QPSK constellation. Fig. [2] shows the 
performance in a 4 x 4 MIMO system with 16-QAM. The performance analysis of a 
MIMO system with different number of antennas employing 8-PSK is shown in Fig. [3l 
In figures [T} [21 and [3l the curved lines with the stars represent the performance of the 
system using relaxation model (|32l) , while a simple rounding algorithm, as described in 
Section |Wl transforms matrix Y to the binary vector u. The ML decoding performance 
is also denoted by a curved line with circles. By increasing the dimension, the resulting 
gap between the relaxation model (l32l) and the ML decoding increases. However, using 
the randomization Algorithm I with M^and = 30 to 50 significantly decreases this gap 
(curved line with diamonds). The curved lines with squares show the performance of the 
relaxation model (|34l) with a simple rounding, in which all the non-negative constraints 
are included. This curve is close to ML performance. It is clear that the relaxation model 
(|34l ) is much stronger than the relaxation model (|32l) . Note that adopting different number 
of non-negative constraints will change the performance of the system between the two 
curves with diamonds and squares. In other words, the trade-off between complexity and 
performance relies on the number of extra non-negative constraints. 

Fig. m compares the two proposed Randomization procedure for the relaxation model 
(|32l ) and (|34l) . The effect of the randomization methods. Algorithm I and II, for the 
relaxation model (|32l) is shown. As expected. Algorithm II performs slightly better. 
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while its computational complexity is lower. The solution of the relaxation model in 
(|34|) . in most cases, corresponds to the optimal solution of the original problem (fTTI) . 
In the other words, because the model in (l34l) is strong enough, there is no need for 
the randomization algorithm. Several compromises for improving the performance can 
be done, e.g. including only some of the non-negative constraints in (l34l) and/or using a 
randomization procedure with a fewer number of iterations. 

In order to reduce the computational complexity of the proposed method, the LLL 
lattice basis reduction is implemented as a pre-processing step for the relaxation model 
(|34|) . Fig. [5] and Fig. [6] show the effect of using the LLL lattice basis reduction in 2 x 2 
and 4x4 multiple antenna systems with 64-QAM and 256-QAM. In a system with 
64-QAM and 256-QAM, the performance of the relaxation model (|34|) is close to the 
ML performance with K = 8 and K = 16, respectively. By using LLL reduction and 
considering L = \og2{K) symbols around the initial point, the performance degradation 
is acceptable, see Fig. |5] and Fig. [6l Note that the resulting gap in the performance is 
small, while the reduction in computational complexity is substantial. 

NxN Multiple Antenna Systems withi QPSK Modulation 

10° 



10"' 



m 10 
m 

10-= 

10-* 

-5 5 10 15 

Fig. 1. Performance of the proposed model i32i and l l34b in a MIMO system with N transmit and A'^ receive antennas 
employing QPSK 



Model (32) - Rounding 
Model (32) - Algorittim I 
Model (34) - Rounding 
ML Decoding 
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NxN Multiple Antenna Sysmtems with 16QAM Modulation 
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Fig. 2. Performance of the proposed model i32l and 04b in a MIMO system with A'^ transmit and A'^ receive antennas 
employing 16-QAM 



NxN Multiple Antenna Systems with 8PSK Modulation 
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Fig. 3. Performance of the proposed model | |32| | and 04b in a MIMO system with N transmit and A'^ receive antennas 
employing 8-PSK 
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4x4 Multiple Antenna System with 16QAM Modulation 
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Fig. 4. Different randomization algorithms in a MIMO system with 4 transmit and 4 receive antennas employing 
16-QAM 




Fig. 5. Performance of using LLL lattice basis reduction for relaxation model l l34b in a 2 x 2 MIMO system with 

L = log2{K) 
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Complexity Reduction Using LLL Lattice Basis Reduction 
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Fig. 6. Performance of using LLL lattice basis reduction for relaxation model l l34b in a 4 x 4 MIMO system with 
L = log2(if) 

B. Complexity Analysis 

Semi-definite programs of reasonable size can be solved in polynomial time within 
any specified accuracy by IPMs. IPMs are iterative algorithms which use a Newton-like 
method to generate search directions to find an approximate solution to the nonlinear 
system. The IPMs converge vary fast and an approximately optimal solution is obtained 
within a polynomial number of iterations. For a survey on IPMs see [47], [48]. In the 
sequel, we provide an analysis for the worst case complexity of solving models (|32l) and 
(|34l) by IPMs. 

It is known (see e.g. [49]) that a SDP with rational data can be solved, within 
a tolerance e, in 0(v^log(l/e)) iterations, where m is the dimension of the matrix 
variable. Note that for the SDP problems and dH]), m = N{K - 1) + 1. 

The computational complexity for one interior-point iteration depends on several 
factors. The main computational task in each iteration is solving a linear system an order 
determined by the number of constraints, p. This task requires 0{p^) operations. The 
remaining computational tasks involved in one interior-point iteration include forming 
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system matrix whose total construction requires 0{pm? + p^m?) arithmetic operations. 
Thus, the complexity per iteration of the IPM for solving SDP problem whose matrix 
variable is of dimension m and number of equality constraints p, is 0{pm? +p^m'^ +p^). 
This means for a given accuracy e, an interior-point method in total requires at most 
0{p{m^ + pm? +p'^)^/m log(l/e)) arithmetic operations. 

Since the SDP relaxation (|32|) contains 0{K'^N) equality constraints, it follows that 
a solution to (|32l) can be found in at most 0(A^^'^i^'^'^ log(l/e)) arithmetic operations. 
SDP relaxation (|34l) contains 0{K'^N) equations and 0{K'^N'^) sign constraints. In order 
to solve relaxation (l34l) . we formulate the SDP model as a standard linear cone program 
(see e.g. [50]) by adding some slack variables. The additional inequality constraints make 
the model in (l34l) considerably stronger than the model in (|32|) (see numerical results), 
but also more difficult to solve. An interior-point method for solving SDP model (|34l) 
within a tolerance e requires at most 0{N^'^ K^'^ \og{l / e)) arithmetic operations. Since 
the problem sizes of interest are moderate, the problem in (l34l) is tractable. However, there 
exist a trade-off between the strength of the bounds and the computational complexity 
for solving these two models (see Section HIH). 

The complexity of the randomization procedure applied to the model (|32l) is negli- 
gible compared to that of solving the problem itself. Namely, if we denote the number of 
randomization iterations by Nrand, then the worst case complexity of the randomization 
procedure is 0{NKNrand)- 

The problems (|32l) and (|34l) are polynomially solvable. These problems have many 
variables; however, they contain sparse low-rank (rank-one) constraint matrices. Exploit- 
ing the structure and sparsity characteristic of semi-definite programs is crucial to reduce 
the complexity. In [51], it is shown that rank-one constraint matrices (similar to our 
problems) reduce the complexity of the interior-point algorithm for positive semi-definite 
programming by a factor of NK. In other words, the complexities of the SDP relaxation 
problems ^ and dH) are decreased to 0{N'^-^K^-^ log(l/e)) and 0{N^-^K^-^ log(l/e)), 
respectively. Also, implementing the rank-one constraint matrices results in a faster 
convergence and a saving in the computation time and memory requirements. It is 
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worth mentioning that when we use the LLL lattice basis reduction, the value of K 
is replaced with L in the aforementioned analysis. As mentioned before, this value is 
much smaller than K, e.g. in our simulation results L = log2(-ft'), which results in 
reducing the computational complexity. 

C. Comparison 

The worst-case complexity of well-known SD method [15], [24] is known to be 
an exponential function of dimension M over all ranges of rate and SNR [24]. The 
complexity analysis shows that our proposed SDP algorithms possess a polynomial-time 
worst case complexity. It should be emphasized that in real time problems, the time spent 
for decoding the received vector is important and it can be considered as a measure of 
the complexity. 

In the following, the worst case complexities of the algorithm based on model 
(|32|) . the method proposed in [30], the method in [27], and the SD algorithm [15] are 
compared with different random values of input vector, channel matrix, and noise for 
Eb/No = {— 5,0,5,10,15}. For each value of E^/Nq, the algorithms are performed for 
10^ times and the maximum time spent for the decoding procedure is saved in MaxTime. 
The average time spent for decoding each case is stored in AveTime. 

It should be emphasized that the MaxTime for each case depends on how the al- 
gorithm is implemented. There are numerous variants for SD algorithm. In the following, 
we have implemented the SD algorithm based on the Schnorr-Euchner strategy proposed 
in [15]. Moreover, the simulations of the proposed algorithms are implemented by one 
of the simplest available packages, the SDPA package [41]. However, by utilizing the 
sparsity of the constraint matrices as suggested in [51] and using the DSDP package, the 
computed AveTime and MaxTime can be reduced dramatically (a factor of NK in the 
analysis), without any performance degradation. 

Table H] shows the simulation results for a MIMO system with M = N = A 
employing 16-QAM. The maximum time for decoding a symbol using SD algorithm 
is much longer than the corresponding time in the proposed SDP relaxation method. The 
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TABLE I 

Comparison of MaxTime for different methods in a MIMO system with 4 transmit and 4 receive 

ANTENNAS EMPLOYING 16-QAM 





-5 





5 


10 


15 


Model m 


0.1037 


0.1095 


0.1108 


0.1178 


0.1196 


Method [30] 


0.0685 


0.0640 


0.0697 


0.0735 


0.0624 


Method [27] 


0.0580 


0.0633 


0.0536 


0.0646 


0.0596 


SD Method 


61.8835 


47.0480 


28.0347 


4.3848 


2.2477 



Other three methods have comparable MaxTime. As it is also shown in the analysis, the 
proposed Model III is more complex compared to the two other SDP methods. However, 
this method outperforms the other SDP methods in [27] and [30]. 

The relaxation model (|32l) outperforms the SDP methods proposed in [27] and 
[29]. Fig. |7] compares the performance of [29] and the relaxation model (|32l ) and the 
performance of the method proposed in [27] is shown in Fig. [8] in a MIMO system with 
4 transmit and receive antennas. The order of the complexity of [27] is comparable to the 
proposed model (l32l) and the order of the complexity of [29] is less than that of the model 
(l32l) (O(iV^) vs. 0{N^-^)). The method in [27] can handle QAM constellations; however, 
it achieves near ML performance only in the case of BPSK and QPSK constellations. 
Also, the method in [29] is limited to PSK constellations. Note that our proposed methods 
can be used for any arbitrary constellation and labeling. 

The comparison of the performance of the relaxation model proposed in [30] and 
that of our method is shown in Fig. [8] (4 x 4 antenna system employing 16-QAM). It 
is observed that the SDP relaxation models (|32|) and (|34|) perform better than [30]. The 
order of the complexity of [30] is the same as that of the model (|24|) . while the model 
dMl) is more complex iO{N^-^) vs. 0{N^-^)). 

Although the worst case complexities of the SD algorithm [15], [24] and the other 
variants are exponential, in several papers, the average complexity of these algorithms are 
investigated. In [26], it is shown that generally, there is an exponential lower bound on 
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4x4 Multiple Antenna Systems with 8PSK Modulation 



* - Model proposed in [29] 
♦ — Model (32) - Algorithm I 
• — Model (34) - Rounding 
■ — ML Decoding 




Fig. 7. Comparison of the relaxation model proposed in [29] and that in our method in a MIMO system with 4 
transmit and receive antennas, employing 8-PSK modulation 



4x4 Multiple Antenna Sysmtems with 16QAM Modulation 




Fig. 8. Comparison of the relaxation model proposed in [30] and that in our method in a MIMO system with 4 
transmit and receive antennas, employing 16-QAM modulation 
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TABLE II 

Comparison of AveTime for different methods in a MIMO system with 4 transmit and 4 receive 

ANTENNAS EMPLOYING 16QAM 



Et/No 


-5 5 10 15 


Model III 
Method [30] 
Method [27] 
SD Method 


0.0372 0.0377 0.0394 0.0428 0.0417 
0.0130 0.0134 0.0142 0.0156 0.0156 
0.0116 0.0118 0.0126 0.0141 0.0141 
0.0449 0.0139 0.0060 0.0026 0.0016 



the average complexity of the SD algorithm. However, it is shown that for large values of 
Eh /No and small values of dimension M, the average complexity can be approximated 
by a polynomial function of dimension M. 

In Table [III the average time AveTime spent for decoding the received vectors in 
the previous scenario is shown. As it can be seen, the average complexity of all SDP 
methods is gradually increasing with E^/Nq while the average complexity of SD method 
is decreasing exponentially. This suggests that for different dimensions M and values 
of Eb/No, there is a threshold that the proposed SDP methods perform better than SD 
algorithm even in terms of the average complexity. However, Table |I] shows that how 
inefficient SD algorithm performs in terms of the worst-case complexity. 

In the following, in Tables UlI] and |Wl the performance of the proposed algorithm 
based on Model III and SD algorithm are shown in terms of AveTime and MaxTime, 
for different number of antennas and constellations. It can be seen that, in terms of the 
worst-case complexity the proposed algorithm based on Model III always outperforms 
SD algorithm. Generally, we can conclude that by increasing the dimension and rate, 
the range of E^/Nq that the proposed model outperforms the SD algorithm increases. In 
order to show that the MaxTime values are not sporadic, we also provide the values 
of AveMaxTime in Table HVl This number is the average of the largest 100 decoding 
times in each case. 

The performance of the proposed SDP relaxation model (l34l) . Model IV, is close to 
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TABLE III 

Decoding Time in a MIMO system with 4 transmit and 4 receive antennas employing QPSK 





Eb/No 


-5 





5 


10 


15 




Model III 


0.0154 


0.0156 


0.0238 


0.0278 


0.0236 


AveTime 
















SD Method 


0.0199 


0.0074 


0.0046 


0.0028 


0.0020 


MaxTime 


Model III 
SD Method 


0.4271 
28.326 


0.4251 
26.3109 


0.4765 
25.4260 


0.7572 
2.2232 


0.8417 
0.9663 



TABLE IV 



Decoding Time in a MIMO system with 8 transmit and 8 receive antennas 





Et/No 


-5 





5 


10 


15 




AveTime 


Model III 
SD Method 


0.0152 
0.6005 


0.0152 
0.1061 


0.0174 
0.0319 


0.0224 
0.0149 


0.0306 
0.0052 


QPSK 


MaxTime 


Model III 
SD Method 


0.0965 
433.3972 


0.0655 
179.0310 


0.1666 
19.7889 


0.6586 
16.7787 


0.6959 
7.7819 




AveMaxTime 


Model III 
SD Method 


0.0658 
73.4830 


0.0587 
16.3274 


0.0642 
6.1652 


0.1492 
5.9249 


0.2109 
1.8074 




AveTime 


Model III 
SD Method 


0.0936 
42.2894 


0.0948 
1.6575 


0.0984 
0.4762 


0.1050 
0.2955 


0.1059 
0.1080 


I6-QAM 


MaxTime 


Model III 
SD Method 


0.2867 
8633.8 


0.2974 
383.40 


0.2772 
290.89 


0.2916 
121.92 


0.3273 
91.032 




AveMaxTime 


Model III 
SD Method 


0.1574 
411.6743 


0.1580 
15.3724 


0.1633 
4.5987 


0.1682 
2.9336 


0.1712 
1.0162 



the ML performance. Similar to the SDP relaxation model (|32l) . the algorithm based 
on Model IV outperforms the SD algorithm in terms of the worst case complexity 
(polynomial vs. exponential). Furthermore, by using the LLL lattice basis reduction before 
the proposed SDP model, the complexity is reduced, with an acceptable degradation in 
the performance (as shown in Simulation Results section). 

As a final note, we must emphasis that in the complexity analysis for model (l34l) . 
we have considered all the non-negative constraints. This suggests that the complexity 
of this model is not tractable. However, it is not required to consider all the non- 
negative constraints. In order to implement this model more efficiently, we can solve 
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the SDP relaxation (|32|) . and solve the SDP relaxation (|34|) with only the most violated 
constraints. These constraints correspond to those positions in matrix Y where their 
values are the minimum negative numbers. Implementing Model IV based on the most 
violated constraints reduces the complexity to almost a number of times more complex 
compared to the Model III. 

VIII. Conclusion 

A method for quasi-maximum likelihood decoding based on two semi-definite re- 
laxation models is introduced. The proposed semi-definite relaxation models provide a 
wealth of trade-off between the complexity and the performance. The strongest model 
provides a near-ML performance with polynomial-time worst-case complexity (unlike the 
SD that has exponential-time complexity). Moreover, the soft decoding method based on 
the proposed models is investigated. By using lattice basis reduction the complexity of 
the decoding methods based on theses models is reduced. 
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Appendix I 
Lagrangian Duality 

In this appendix, we show that Lagrangian duality can be used to derive the SDP 
relaxation problem (|24l) . We first dualize the constraints of (fTTI) . and then derive the SDP 
relaxation from the dual of the homogenized Lagrangian dual. Finally, we project the 
obtained relaxation onto the minimal face. The resulting relaxation is equivalent to the 
relaxation (l24l) . 

Consider the minimization problem in (fTTI) . According to [31], for an accurate 
semi-definite solution, zero-one constraints should be formulated as quadratic constraints. 
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Therefore, 

mill u^Qu + 2c^u 
s.t. ||Au - ejvlP = 

uf — Ui \/i — 1, ■ ■ ■ ,n. (43) 

First, the constraints are added to the objective function using lagrange multipliers A and 

w = [wi, ■ ■ ■ ,Wnf: 

fio = minmax ju^Qu + 2c^u 

u A,w 

+ A (u^A^Au - 2e^Au + e^ejv) 

n 

+ ^Wi{u]-Ui)}. (44) 
Interchanging min and max yields 

pf'O ^ pi-c — max min ju^Qu + 2c^u 

A,w u 

+ A (u^A^Au - 2e^Au + e^e^v) 

n 

1=1 

Next, we homogenize the objective function by multiplying t with a constrained scalar 
uo and then increasing the dimension of the problem by 1. Homogenization simplifies 
the transition to a semi-definite programming problem. Therefore, we have 

l^o > l^c = max min {u'' [Q + AA^A + Diag(w)] u 

A,w u,ug=l 

- (2Ae^A - 2c^ + w^) uqU 

+ XeJfeN } , (46) 

where Diag(w) is a diagonal matrix with w as its diagonal elements. Introducing a 
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Lagrange multiplier wq for the constraint on uq, we obtain the lower bound fin 



/^o > > yU7e=max min {u^TQ + AA^A + Diag(w)lu 

A,w,uio uo,u 

- {2XefjA - 2c^ + w^) uqu 

+ XeljeNul + Wo {ul - l) } • 



(47) 



Note that both inequalities can be strict, i.e. there can be duality gaps in each of the 
Lagrangian relaxations. Also, the multiplication of Xej^ejy by Uq is a multiplication by 



1. Now, by grouping the quadratic, linear, and constant terms together and defining 
[uo, u^]^ and w"^ = [wq, w"^]^, we obtain 



U 



/i7^= max min {u"^[£q + Arrow (w) + A£a]u — Wo] , 



(48) 



where 



-A^eAT A^A 



N 



-ex ® Bat Iat (g) {exef^) 



Arrow(w) 
and £q = 



Wo 

1 



2Wl:„ 
T 



1 T 

Diag(wi,„) 



(49) 





c Q 

Note that we will refer to the additional row and column generated by the homog- 
enization of the problem as the 0-th row and column. There is a hidden semi-definite 
constraint in (l48l) . i.e. the inner minimization problem is bounded below only if the 
Hessian of the quadratic form is positive semi-definite. In this case, the quadratic form 
has minimum value 0. This yields the following SDP problem: 

max — Wo 

s.t. Cq + Arrow(w) + A^a h 0. (50) 

We now obtain our desired SDP relaxation of (|43l) as the Lagrangian dual of (l50l) . By 
Introducing the (n + 1) x (n + 1) dual matrix variable Y ^ 0, the dual program to the 



SUBMITTED TO IEEE TRANS. ON INFO. THEORY 



35 



SDP (gOl) would be 



mill trace CqY 
s.t. diag(Y) = (l,Yo,i:„)^ 



traceY/^A 







Y ^ 



(51) 



where the first constraint represents the zero-one constraints in (l43l) by guaranteeing that 
the diagonal and 0-th column (row) are identical (matrix Y is indexed from 0); and the 
constraint Au = e^r is represented by the constraint trace Y£a = 0. Note that if the 
matrix Y is restricted to be rank-one in (|5T1) . i.e. 



for some u G M", then the optimal solution of (fSTI) provides the optimal solution, u, for 



Since the matrix Cx ^ is a positive semi-definite matrix; therefore, to satisfy the 
constraint in (|5T1) . Y has to be singular. This means the feasible set of the primal problem 
in (ISTI) has no interior [34] and an interior-point method may never converge. However, 
a simple structured matrix can be found in the relative interior of the feasible set in order 
to project (and regularize) the problem into a smaller dimension. 

As mentioned before, the rank-one matrices are the extreme points of the feasible 
set of the problem in (|5TI) and the minimal face of the feasible set that contains all these 
points shall be found [34]. 

From Theorems [2] and [3l we conclude that Y ^ is in the minimal face if and 
only if Y = VRV^, for some R ^ 0. By substituting VRV^ for Y in the SDP 
relaxation (|5TI) . we get the following projected SDP relaxation which is the same as the 



1 




Y 



u 



(El). 
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SDP relaxation in (|24l) : 

fiRi = min trace (V^£qV)R 



s.t. diag(VRV ) = (1, (VRV ;o,i:nj 
R ^ 0. 



(52) 



Note that the constraint trace(V-'"£AV)R = can be dropped since it is always satisfied, 
i.e. CxV = 0. 



Appendix II 
Proofs 

A. Lemma [7] 

Let X G Mkxn and e^^^X = e^. Since ^kx(k~i) is a Kx. {K—1) matrix containing 
a basis of the orthogonal complement of the vector of all ones, i.e., V^^^^^^-je/^ = 0, 
and 

e^Fi^xAf = e^, (53) 

we have 



X — Fa'xAT + Vii'x(X-l)Z, 



where Z G M.{K-i)>iN- From 



O(i^-i) 



and 



Vi<:x(x-i)Z 



(54) 



(55) 



(56) 



it follows that Z = X(l : - 1), 1 : A^). 
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B. Theorem |2] 

Let Y e ^ be an extreme point of JF, i.e. 

1 



XX 



(57) 



for some x = vec(X), X G Srxn- From Lemma [H it follows that every matrix X G 
Skxn is of the form X = Fkxn + ^kx{k-i)^ where X = X(l : K — 1,1 : N). From 
the properties of the Kronecker product (see [32]), it follows 



X = vec(X) = — (eA'AT 



fl 



N 



^Kx{K-1))G{K-1)n) 



N 



V 



where x = vec(X). Let p-^ 
1 



1 i 



Kx{K-l))^, 



and 



W := 



K 



GkN — (ItV ® ^Kx(K~1))G(K-1)xn), lAf ® Vxx(A'-l) 



(58) 



(59) 



Therefore, x = Wp, and 





1 


p^W^ 


Y = 






Wp 


Wpp^W^ 



VRV^, 



(60) 



where R := pp^, i.e. 



R 



x-* 



rp 

XX 



y 0. 



(61) 



Since x is a binary vector, it follows that r^j G {0,1}, Vi,j G {0, . . . , N{K — 1)}, 
and diag(xx^) = x. The proof follows analogously for any convex combination of the 
extreme points from JF. 



C. Theorem 13 

Fix u G W and let 



1 
u 
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Considering the constraint on this vector in (fT4l) . we divide u into sub-vectors of 
length K. In each sub-vector all the elements are zero except one of the elements which 
is one. Therefore, there are K'^ different binary vectors in the set U. Consider the entries 
of the 0-th row of Y. Note that y^j = 1 means that the j-th element of u is 1. In addition, 
there is only one element equal to 1 in each sub-vector. Therefore, there are K^^^ such 
vectors, and the components of the 0-th row of Y are given by 



K 



N-l 



1 



Now consider the entries of Y in the other rows, K . . 



1) If i = j, then, ?/j , = 1 means that the i-th element of the vector u is 1 and there 



are ^ such vectors; therefore, the diagonal elements are 



yi,i 



K 



N-l 



1 

li' 



2) \fi = K{p-l)+q, J = Kip-l)+r,q=^r,q,r e {I, ■ ■ ■ , K},p E {1, ■ ■ ■ ,N}, 
i.e. the element is an off-diagonal element in a diagonal block, then, yij = 1 means 
that the z-th and the j-th elements of u in a sub-vector should be 1 and this is not 
possible. Therefore, this element is always zero. 



3) Otherwise, we consider the elements of the off-diagonal blocks of Y. Then, y, 



1 



means that the i-th and the j-th elements of u in two different sub-vectors are 1 
and there are K'^^^ such vectors; therefore, the elements of the off-diagonal blocks 
are 



1 



K 



N-2 



1 



This proves the representation of Y in (i) in Theorem [3l 
It can be easily shown that 



Y 



1 


of 




In 



1 


1 T 
—P 




On 


In 



1 


of 


On 


w 



(62) 



where W = -^In ® {KIk - ^k)- Note that rank(Y) = 1 + rank(W). 

The eigenvalues of Ia? are 1, with multiplicity and the eigenvalues of KIk — Ea' are 

K, with multiplicity K — 1, and 0. Note that the eigenvalues of a Kronecker product are 
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given by the Kronecker product of the eigenvalues [32]. Therefore, the eigenvalues of 
W are with multiplicity N{K — 1), and 0, with multiplicity A^. Therefore, we have 

rank(Y) = 1 + rank(W) = N{K - 1) + 1. (63) 

This proves (ii) in Theorem [3l 

By (|62|) and (|63l) . we can easily see that the eigenvalues of Y are — , with mul- 

K + N 

tiplicity N{K — 1), — — — , and 0, with multiplicity A^. This proves (iii) in Theorem 

K 

m 

The only constraint that defines the minimal face is Au = b. By multiplying of 
both sides by and using the fact that u is a binary vector, we obtain 

Tuu^ = (diag(uu^))^ . (64) 

This condition is equivalent to 

TY = 0. (65) 

Note that rank(T) = A^. Therefore, we have 

Ar(Y) = {u:ue 7^(T^)}. 

This proves (iv) in Theorem [3l 

Since rank(V) = N{K — 1) + 1 and using Theorem |2l the columns of V span the 
range space of Y. This proves (v) in Theorem [3l 

Appendix III 
Randomization Algorithm 

In this appendix, we present two randomization algorithms to transform Y to a 
binary rank-one matrix. 



SUBMITTED TO IEEE TRANS. ON INFO. THEORY 



40 



A. Algorithm I 

Goemans and Williamson [52] introduced an algorithm that randomly transforms 
an SDP relaxation solution to a rank-one solution. This approach is used in [27] for the 
quasi maximum likelihood decoding of a PSK signalling. This technique is based on 
expressing the BPSK symbols by {—1,1} elements. After solving the relaxation problem 
in [27], the Cholesky factorization is applied to the n x n matrix Y and the Cholesky 
factor V = [vi, . . . , v„] is computed, i.e. Y = VV"^. In [27], it is observed that one 
can approximate the solution of the distance minimization problem, u, using V, i.e. Ui 
is approximated using Vj. Thus, the assignment of —1 or 1 to the vectors {vi, . . . , v„} 
is equivalent to specifying the elements of u. 




Fig. 9. Representation of the randomization algorithm in [52] 

It is shown that norms of the vectors {vi, . . . , v„} are one, and they are inside 
an n-dimensional unit sphere [27], see Fig. |9l These vectors should be classified in 
two different groups corresponding to 1 and —1. In order to assign —1 or 1 to these 
vectors, the randomization procedure generates a random vector uniformly distributed 
in the sphere. This vector defines a plane crossing the origin. Among given vectors Vj, 
i = 1, . . .n, all the vectors at one side of the plane are assigned to 1 and the rest are 
assigned to —1, as shown in Fig. |9l This procedure is repeated several times and the 
vector u resulting in the lowest objective function is selected as the answer. 

In our proposed approach, the variables are binary numbers. In order to implement 
the randomization procedure of [52], we bijectively map the computed solution of the 
{0, 1} SDP formulation to the solution of the corresponding { — 1, 1} SDP formulation. 
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More precisely, we use the following mapping: 





1 


■■ 


• " 




-1 


2 ■■ 


■ 


M = 










-1 


■■ 


■ 2 



Y{_i,i| = MY{o,i}M^, (66) 

where Y{o,i} is the resulting matrix from the relaxation model (|32l) or (|34|) and Y{_i i} is 
its corresponding matrix with { — 1, 1} elements. Using (l66l) . the solution for (fTT)) can be 
computed using a similar randomization method as in [27]. The computational complexity 
of this randomization algorithm is polynomial [27]. 

Considering zero-one elements in our problem, we propose a new randomization 
procedure inspired by [52]. This algorithm can be applied to {0, 1} formulation directly. 
Therefore, the complexity of the whole randomization procedure is reduced, since the 
preprocessing step, i.e. bijective mapping in (|66l) . is omitted. 

B. Algorithm II 

After solving the relaxation model (|32|) or (|34|) . the Cholesky factorization of Y 
results in a matrix V = [vi, . . . , Vn] such that Y = VV^. The matrix Y is neither 
binary nor rank-one. Therefore, norms of the resulting vectors Vj are between zero and 
one. These vectors are depicted in Fig. \T0\ Intuitively, a sphere with a random radius 
uniformly distributed between zero and one has the same functionality as the random 
plane in Fig. |9l 

In order to assign or 1 to these vectors, the randomization procedure generates 
a random number, uniformly distributed between and 1, as the radius of the sphere. 
Among given vectors Vj, i = l,...,n, all the vectors whose norms are larger than 
this number are assigned to 1 and the rest are assigned to 0. In another variation of this 
algorithm, the radius of the sphere can be fixed, and norms of these vectors are multiplied 
by a random number. This procedure is repeated several times and the vector u resulting in 
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Fig. 10. Graphic representation for the proposed randomization algorithm 

the smallest objective function value in (fTTl) is selected as the solution. Simulation results 
confirm that the proposed method results in a slightly better performance for the lattice 
decoding problem compared to the first algorithm. Also, the computational complexity 
of the randomization algorithm is decreased, due to the removal of the preprocessing 
step in (|66l) . It is worth mentioning that, according to (1571 ) and (|6T| ). the randomization 
procedure can be implemented for the matrix R, which results in further reduction in the 
computational complexity. 
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